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Abstract. We investigate, by means of extensive Monte Carlo simulations, the magnetic critical behavior of 
the three-dimensional bimodal random-field Ising model at the strong disorder regime. We present results in 
favor of the two-exponent scaling scenario, f) = 2rj, where r\ and fj are the critical exponents describing the 
power-law decay of the connected and disconnected correlation functions and we illustrate, using various 
finite-size measures and properly defined noise to signal ratios, the strong violation of self-averaging of the 
model in the ordered phase. 
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1 Introduction 

Phase transitions in pure systems [I] are already well un- 
derstood. The critical behavior of all physical quantities 
can be described via critical exponents which are related 
through (hyper) scaling relations to each other, so that 
only two independent exponents remain. In contrast, phase 
transitions in systems with (quenched) disorder [2] exhibit 
many puzzles and are still far from being understood. In 
statistical physics, and after the pioneering work of Imry 
and Ma [3], the random-field Ising model (RFIM) is a 
widely studied prototypical disordered system. It is be- 
lieved [415] to be in the same universality class as the 
diluted antifcrromagnet in a field, which can be studied 
experimentally [BJ. 

Originally, it was believed [71819] that the critical be- 
havior of the d-dimensional RFIM is equal to that of the 
d — 2 pure ferromagnet, implying that the d = 3 RFIM 
exhibits no ordered phase. This is of course not true, as 
has been shown rigorously [lOj a few years later In 
the meantime, a scaling theory [12113114] for the RFIM 
was developed, where the dimension d has been replaced 
by d — 8 in the hyperscaling relations, 8 being a third in- 
dependent critical exponent, in contrast to the pure case. 
An alternative approach |15ll6j leads to the consequence 
that 8 is not independent but related to the exponent r\ 
via 8 = 2 — rj (two-exponent scaling scenario). Further evi- 
dence for the existence of only two independent exponents 
was found by high-temperature expansions [17j . This was 
mainly confirmed for the RFIM with a Gaussian random- 
field distribution in d — 3 by Monte Carlo (MC) simula- 
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tions [18] and exact ground-state calculations [19 20121] . 
and in d = 4 also by ground-state calculations [22] , The 
first goal of the present paper is the numerical verifica- 
tion of the above scenario for the case of the RFIM with 
a bimodal random-field distribution, yet having in mind 
that the precision of numerical simulations may never be 
good enough to actually prove that this scenario is truly 
verified; rather we should see this approach as a further 
test in favor of the proposed scaling picture of the RFIM. 
A first attempt towards this direction has been presented 
by Rieger and Young [53] in an early paper, yet their anal- 
ysis was restricted to rather small lattice sizes (L < 16, L 
being the linear lattice dimension). 

The second part of our study focuses on another im- 
portant relevant issue in the study of disordered systems, 
that of self-averaging. Although it has been known for 
many years now that for (spin and regular) glasses there 
is no self-averaging in the ordered phase [24] . for random 
fcrromagnets such a behavior was first observed for the 
RFIM in a paper by Dayan et al. [25] and some years later 
for the random versions of the Ising and Ashkin- Teller 
models by Wiseman and Domany [2BJ. These latter au- 
thors suggested a finite-size scaling (FSS) ansatz describ- 
ing the absence of self-averaging and the universal fluc- 
tuations of random systems near critical points that was 
refined and put on a more rigorous basis by the intriguing 
renormalization-group work of Aharony and Harris [27] , 
Ever since, the subject of breakdown of self-averaging is an 
important aspect in several theoretical and numerical in- 
vestigations of disordered spin systems [28 29 30 31 32 33 34 35 36137 
In fact, most recently, Efrat and Schwartz [39] showed that 
the property of lack of self-averaging in disordered systems 
may be turned into a useful tool that can provide an inde- 
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pendent measure to distinguish between the ordered and 
disordered phases of the system. In view of this increasing 
interest in the understanding of the self-averaging prop- 
erties of disordered systems we attempt in the present 
paper to examine and apply the theoretical predictions to 
the RFIM, a suitable candidate which is already known 
to suffer from a strong violation of self-averaging in terms 
of the thermal properties [33136] . 

The Hamiltonian describing the RFIM is 

H = -JJ2 SiSj - hJ^hiSi, (1) 

<i,j> i 

where Si are Ising spins, J > is the nearest-neighbors 
ferromagnetic interaction, h is the disorder strength, also 
called randomness of the system, and hi are independent 
quenched random fields obtained here from a bimodal dis- 
tribution of the form 

P{hi) = \[5{hi-l) + 8(hi + l]. (2) 

The critical value of h, denoted as critical disorder strength 
h c , above which no phase transition occurs in the ther- 
modynamic limit is known with good accuracy for the 
bimodal RFIM to be h c = 2.21(1) [19I40J . In particu- 
lar, in the present paper we consider the model defined 
in equations ^ and ([2]) in the strong disorder regime 
(h 6 [1.9,2.25]) and embedded in simple cubic lattices 
with periodic boundary conditions and N — L 3 spins, 
where the linear dimension L takes values in the range 
L = 4 — 32. The ensembles of random-field realizations 
Q (q — 1,2, ... ,Q) simulated for the present study are 
the following: Q L < 16 = 1000 and Q L >xe = 500. The total 
computer time used for the simulations was several (~ 8) 
Intel Pentium-Pro months. 

The rest of the paper is laid out as follows: In Sec- 
tion [2] we briefly describe here our numerical implementa- 
tion, that is based on an efficient implementation of the 
Wang-Landau (WL) method [H] , as also presented in a 
number of papers dealing with the simulation aspects of 
disordered systems [42143144] . In Section [3] we discuss the 
FSS behavior of the magnetic quantities of the model and 
provide evidence in favor of the two-exponent scaling sce- 
nario. Subsequently, in Section @] we investigate the self- 
averaging properties of the model in terms of broad distri- 
butions and we illustrate using various finite-size measures 
and properly defined noise to signal ratios, the strong vio- 
lation of self-averaging of the model in the ordered phase. 
Finally, we summarize our conclusions in Section [5] 



2 Numerical approach 

MC simulations in statistical physics now have a history 
of nearly half a century starting with the seminal work of 
Metropolis. While the Metropolis algorithm has been es- 
tablished as the standard algorithm for importance sam- 
pling it suffers from two problems: the inability to directly 
calculate the partition function, free energy or entropy, 



and critical slowing down near phase transitions and in 
disordered systems. In a standard MC algorithm a series of 
configurations is generated according to a given distribu- 
tion, usually the Boltzmann distribution in classical sim- 
ulations. While this allows the calculation of thermal av- 
erages, it does not give the partition function, nor the free 
energy. They can only be obtained with limited accuracy 
as a temperature integral of the specific heat, or by using 
maximum entropy methods. The problem of critical slow- 
ing down has been overcome for second-order phase transi- 
tions by cluster update schemes [45) . For first-order phase 
transitions and systems with rough free-energy landscapes 
(such the present model under study) a decisive improve- 
ment was achieved some years ago via the WL algorithm [?T] 

In the last few years our group has used an entropic 
sampling implementation of the WL algorithm [JT] to sys- 
tematically study some simple [42!, but also some more 
complex systems [36 40 43 44] . One basic ingredient of this 
implementation is a suitable restriction of the energy sub- 
space for the implementation of the WL algorithm. This 
was originally termed as the critical minimum energy sub- 
space (CrMES) restriction [52] and it can be carried out 
in many alternative ways, the simplest being that of ob- 
serving the finite-size behavior of the tails of the energy 
probability density function (e-pdf) of the system [4"2"] . 

Complications that may arise in complex systems, i.e. 
random systems or systems showing a first-order phase 
transition, can be easily accounted for by various simple 
modifications that take into account possible oscillations 
in the e-pdf and expected sample-to-sample fluctuations 
of individual realizations. In our recent papers [40143144] . 
we have presented details of various sophisticated routes 
for the identification of the appropriate energy subspace 
(Ei,E2) for the entropic sampling of each realization. In 
estimating the appropriate subspace from a chosen pseu- 
docritical temperature one should be careful to account 
for the shift behavior of other important pseudocritical 
temperatures and extend the subspace appropriately from 
both low- and high-energy sides in order to achieve an ac- 
curate estimation of all finite-size anomalies. Of course, 
taking the union of the corresponding subspaces, insures 
accuracy for the temperature region of all studied pseud- 
ocritical temperatures. 

The up to date version of our implementation uses a 
combination of several stages of the WL process. First, we 
carry out a starting (or preliminary) multi-range (multi- 
R) stage, in a very wide energy subspace. This prelimi- 
nary stage is performed up to a certain level of the WL 
random walk. The WL refinement is G(E) — > / • G(E), 
where G{E) is the density of states (DOS) and we follow 
the usual modification factor adjustment ,fj+i = \[Tj and 
fx = e. The preliminary stage may consist of the levels 
: j = 1, . . . ,j = 18 and to improve accuracy the process 
may be repeated several times. However, in repeating the 
preliminary process and in order to be efficient, we use 
only the levels j — 13, . . . , 18 after the first attempt, us- 
ing as starting DOS the one obtained in the first random 
walk at the level j = 12. From our experience, this prac- 
tice is almost equivalent to simulating the same number of 
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independent WL random walks. Also in our recent studies 
we have found out that is much more efficient and accurate 
to loosen up the originally applied very strict flatness cri- 
teria [H]. Thus, a variable flatness process starting at the 
first levels with a very loose flatness criteria and assuming 
at the level j = 18 the original strict flatness criteria is now 
days used. After the above described preliminary multi-R 
stage, in the wide energy subspace, one can proceed in 
a safe identification of the appropriate energy subspace 
using one or more alternatives outlined in reference [42] . 

The process continues in two further stages (two-stage 
process), using now mainly high iteration levels, where 
the modification factor is very close to unity and there is 
not any significant violation of the detailed balance con- 
dition during the WL process. These two stages are suit- 
able for the accumulation of histogram data (for instance 
energy- magnetization histograms), which can be used for 
an accurate entropic calculation of non-thermal thermo- 
dynamic parameters, such as the order parameter and its 
susceptibility [42 . In the first (high-level) stage, we follow 
again a repeated several times (typically ~ 5 — 10) multi- 
R WL approach, carried out now only in the restricted 
energy subspace. The WL levels may be now chosen as 
j = 18, 19, 20 and as an appropriate starting DOS for the 
corresponding starting level the average DOS of the pre- 
liminary stage at the starting level may be used. Finally, 
the second (high-level) stage is applied in the refinement 
WL levels j — ji, ■ ■ ■ ,ji + 3 (typically ji — 21), where we 
usually test both an one-range (one-R) or a multi-R ap- 
proach with large energy intervals. In the case of the one-R 
approach we have found very convenient and in most cases 
more accurate to follow the Belardinelli and Pereyra [M] 
adjustment of the WL modification factor according to the 
rule In / ~ i -1 , where t refers to the MC time of the pro- 
cess. Finally, it should be also noted that by applying in 
our scheme a separate accumulation of histogram data in 
the starting multi-R stage (in the wide energy subspace) 
offers the opportunity to inspect the behavior of all ba- 
sic thermodynamic functions in an also wide temperature 
range and not only in the neighborhood of the finite-size 
anomalies. 

Closing this review on the numerical part of this work, 
we would like to point out that, the above scheme may also 
be extended to accumulate multi-parametric histograms, 
and provide information for the critical properties of the 
system that are not easily accessible with other method- 
ologies. For more details on this extension we refer the 
reader to reference [30] . 



3 Two-exponent scaling 

We start by studying the critical behavior of the suscep- 
tibility x an d disconnected susceptibility Xdis, which are 
given, in general, by the standard expressions 

1 N 

X= ¥f J2((S i S j )-(S i )(S j )), (3) 



and 

1 N 

Xdu = J f^,(S i )(S j ), (4) 

where (...) denotes a thermal average. We now follow the 
scaling arguments based on the droplet theory [12113114] . 
as presented by Dayan et al. [25] . According to these argu- 
ments, most samples in a large ensemble of realizations of 
the random field have only one thermally populated min- 
imum. For these samples (called hereafter as typical) the 
two terms in equation ([3]) almost cancel at T c , and while 
each term has separately the divergence of Xdis, he. 

Xdls ~L^ { A- n = *,/v), (5) 

their difference is given by 

Xt ~L 2 -" (2-17 = 7/1/), (6) 

with probability p t ~ 1 (the subscript t stands for typi- 
cal). The cancellation of the most divergent terms in the 
expression of the susceptibility is rigorously true on av- 
erage, according to the Schwartz-Soffer inequality [15] , 
Thus, following reference [15] . one shows the inequality 
4 — 77 < d — 2/3/ 1/, and assuming that this holds as an 
equality (see also reference [55]) one derives the modified 
hyperscaling relation 

4-fj = d-2P/v, (7) 

which can be also written in the form 

(d- 9)v = 2/3 + 7- (8) 

The exponent 9 in the above equation accounts for the vi- 
olation of the hyperscaling relation and is given, according 
to the above equations (Q and (J5J), by 

9 = 2 - 77 + ?y. (9) 

Furthermore, rare samples with more than one mini- 
mum (thermally populated) may appear with a probabil- 
ity p r ~ L~ 6 (where the subscript r stands for rare) and, 
for these samples, the two terms in equation §5§ do not 
cancel at T c , and their difference is of the same order of 
each term, i.e. 

Xr ~ L^. (10) 

Although with increasing L the probability of obtaining 
one of these rare samples decreases, the susceptibility of a 
rare sample gets much larger, relatively to the value of a 
typical sample (xr ^ Xt, see Figure [lj. Thus, we need to 
average over large ensembles of random fields in order to 
obtain good statistics. 

Averaging over Q realizations of the random field and 
using equations ([OJ and (|T0j) , we see that for the first mo- 
ment [x]av (where [. . ] av denotes randomness averaging), 
both typical and rare samples give a comparable contri- 
bution, so that at T c 

[xU = E ~ 1 • ^ + L- e ■ L 4 -' 1 ~ L 2 -'\ (11) 

i=t,r 
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Fig. 1. Magnetic susceptibility as a function of temperature of 
a typical (solid line) and a rare (dotted line) sample for h = 2 
and L = 24. 



The behavior of the average disconnected susceptibility at 
T c is, according to equation ([5]), 



dis\ av 



(12) 



At this point we shall assume that the above power 
laws are valid also for the pseudocritical temperatures of 
some relevant thermodynamic quantities. In particular, we 
will consider here (mainly) the sequence of pseudocritical 
temperatures that correspond to the maxima (for different 
L) of the average susceptibility ([x]™) defined as 



Q 



[xYo 



max{[x] tj 



max 

T 



(13) 



where the subscript q refers to a particular random real- 
ization of the quenched disorder and \q to the correspond- 
ing magnetic susceptibility. These pseudocritical tempera- 
tures will be denoted in the sequel as lij^. (for the case 
of the specific heat the corresponding pseudocritical tem- 
peratures will be denoted as T = T L ^ C ]' )■ In Figure [U[a) 
we plot the maxima of the average susceptibility versus 
the lattice size L in a log- log plot for the case h = 2. The 
full line is a least-squares straight-line fitting which gives 
a slope of 1.49(2), that is according to equation (fTTj) an 
estimate for 2 — rj, yielding 



X] = 0.51(2). 



(14) 



The average disconnected susceptibility at T L ^ X ]* , de- 
noted as [xdis]a«j is depicted in Figuredjb), also for h = 2, 
and is expected to obey the power law (fl~2"j) . Using again 
a least-squares fitting of the data to a straight line in the 
log-log plot we get the estimate for 4 — 77 to be 2.99(3), 
and therefore 

77 = 1.01(3). (15) 



From equations ([14]) and ([15]) we see that the Schwartz- 
Soffcr inequality fj < 2r] [TS] is fulfilled within error bars as 
an equality, supporting the two-exponent scaling scenario, 



h = 2;L = 8-32 
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Fig. 2. FSS in a double- logarithmic scale of (a) the maxima 
of the average susceptibility and (b) the average disconnected 
susceptibility at the temperature where the average suscepti- 
bility attains its maximum, for h = 2 and L — 8 — 32. In both 
panels the solid lines show linear fitting attempts. 



i.e. rj = 2rj [17115116] . The exponent 9 is estimated from 
equation © to be 



9 = 1.50(3). 



(16) 



Overall, the estimated values for rj, r], and 9 are in very 
close agreement with the values 



77 = 0.5 77 = 1 61 = 1.5, 



(17) 



that can be predicted from elementary considerations in 
the ordered phase (see also the discussion in reference [2"5] ) 
and agree with the early dimensional reduction prediction 
d —> d/2 of Shapir [37]. Moreover, the resulting value for 
the exponent /? is very close to zero, as already noted by 
Dayan et al. [IS] . The resulting very small value for j3 (j3 = 
0.007(5)) with a relatively large error could be compared 
with renormalization group results (f3 = 0.0200(5)) [48] 
and zero temperature studies ((3 = 0.017(5)) [2Tj . 

Up to this point, our analysis and discussion followed 
the assumption of a continuous phase transition, which is 
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actually the most widely accepted scenario in the random- 
field community. However, the estimated value for the crit- 
ical exponent (3 raises some doubts of whether this is ac- 
tually true, especially for the present case of the bimodal 
RFIM, for which mean-field theory predicted a tricritical 
point at high values of the random field [35]. This main 
issue has regained interest after the recent observations, 
in both the bimodal [50151) and Gaussian [37152) cases, of 
first-order-like features at the strong disorder regime. In 
particular first-order-like features, such as the appearance 
of the characteristic double-peak structure of the canoni- 
cal e-pdf, have been recently reported for both the Gaus- 
sian and the bimodal distributions of the d — 3 RFIM. 

However, in a recent paper Fytas et al. [43 presented 
a detailed study of the first-order transition features - in 
terms of the surface tension and latent heat - of the d = 3 
RFIM, having in mind that a mere observation of a first- 
order structure is not sufficient for the identification of 
the transition. As it was shown in reference [33], this is 
especially true for the RFIM, since its critical behavior is 
obscured by strong and complex finite-size effects, involv- 
ing also the important issue of the lack of self-averaging, 
discussed here in the next Section. The results of refer- 
ence 43 for two values of the disorder strength in the 
strong disorder regime (h — 2 and h — 2.25) clearly indi- 
cated that the interface tension vanishes and the two peaks 
of the e-pdf move together in the thermodynamic limit 
and therefore provided convincing evidence that the tran- 
sition is continuous and that there in no tricritical point 
along the phase transition line. These results pointed to an 
unconventional continuous transition, in which the e-pdf 
approaches two delta functions that move together in the 
thermodynamic limit. Such an unconventional behavior 
has been originally predicted by Eichhorn and Binder [55] , 
for the the order-parameter pdf of the d = 3 random-field 
three-state Potts model and found its theoretical justifi- 
cation in the framework of a refined FSS theory a few 
months ago in a profound paper by Vink et al. [53] . 

Closing this Section, we may also calculate the critical 
exponents 7 of the susceptibility and 7 of the disconnected 
susceptibility 

7 = 1/(2-17) = 1.95(27) 
7 = v{A-ff) = 3.92(54), (18) 

which also compare favorably with previous estimations 25 2$ 17 W 20'^. a | ue ' \- e ^ 
Note that in equation ([IS]) the large error bounds come 
mainly from the error bounds in v, since we have used the 
value v = 1.31(18) estimated in our previous work [36] . 
We would like to note here that, an analogous study of 
X and Xdis was performed by Eichhorn and Binder [28) 
for the d — 3 three-state random-field Potts model. Their 
findings provided also qualitative evidence for the two- 
exponent scaling description of random-field systems. 
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Fig. 3. The histograms for the probability distribution P(Xq) 
of the susceptibility maxima for L = 24 (and L = 14 in the 
inset). 



our numerical studies of disordered systems are carried out 
near their critical points using finite samples; each sam- 
ple q is a particular random realization of the quenched 
disorder. A measurement of a thermodynamic property w 
yields a different value for the exact thermal average uu q 
of every sample q. In an ensemble of disordered samples of 
linear size L the values of w q are distributed according to a 
probability distribution P(w q ). The behavior of this distri- 
bution is directly related to the issue of self-averaging. In 
particular, by studying the behavior of the width of P{w q ) 
with increasing the system size L, one may address qual- 
itatively the issue of self-averaging, as has already been 
stressed by previous authors [30] . In general, we charac- 
terize the distribution P{w q ) by its average [w) av and also 
by the relative variance 



V' 



[W 2 ]av 



W ; 



(19) 



4 Broad distributions and lack of 
self-averaging 

Let us start this Section by recapitulating some basic defi- 
nitions and concepts on self-averaging. As discussed above, 



Suppose now that w is a singular density of an extensive 
thermodynamic property, such as M or Xi or the singular 
part of e(= E/N) and C. The system is said to exhibit 
self- averaging if R w ~ > as L — > 00. If R w tends to a 

const ^ as L — > 00, then the 
system exhibits lack of self-averaging. 

The importance of the above concepts has been illus- 
trated by Aharony and Harris [27] and their main conclu- 
sions are summarized as follows: (i) Outside the critical 
temperature: R w = 0. In a finite geometry, the correla- 
tion length £ is finite for T ^ T c and it can be found, 
using general statistical arguments, originally introduced 
by Brout [53], that R w cx (£/£) d -> 0, as L -> 00. This 
is called strong self-averaging, (ii) At the critical tem- 
perature there exist two possible scenarios: (a) models 
in which according to the Harris criterion [55] the dis- 
order is relevant (a p > 0): R w ^ 0. Then, the system at 
the critical point is not self-averaging and (b) models in 
which according to the Harris criterion disorder is irrele- 
vant (a p < 0): R w = 0. In this case R w scales as L a / V , 
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h = 2 ; L= 16 - 32 



Fig. 4. The h = 2 histograms of the probability distribution 
P(T* iq ) for L = 24 (main panel) and L = 32 (right inset). The 
left inset shows the FSS of the sample-to-sample variance of 
the average [T~]<,«. 



where a and v are the critical exponents of the pure sys- 
tem, which are the same in the disordered one. This is 
called weak self-averaging, (iii) The pseudocritical tem- 
peratures of the disordered system are distributed 
with a width 6([T^] av ) which scales with the system size as 
5([T^] av ) ~ L~ n , where n = d/2 or n — 1/v, depending on 
whether the disordered system is controlled by the pure or 
the random fixed point, respectively. The above behavior 
is now well established by the pioneering works of Aharony 
and Harris [27] and Wiseman and Domany 26 3D]. 

We start the presentation of our results with the prob- 
ability distribution of the susceptibility maxima P(x* q )- 
P{x* q ) is expected to be a broad distribution with a most 
probable value of the order L 2 ~'' and a tail extending to 
much larger values of the order L 4- '' [55]. Figure [3] and 
the corresponding inset show the normalized histograms 
of P(x* q ) at ft = 2 for lattice sizes L — 24 and L = 14, 
respectively. In both cases P{x* q ) nas its peak at a rather 
small value of %* but there is a clearly developing long 
tail extending to much larger values of %*. As expected, 
the tail is longer for the larger lattice size. Subsequently, 
in the main panel and right inset of Figure 2] we plot the 
h = 2 histograms for the probability distribution P(T* ) 
of the pseudocritical temperatures of the magnetic suscep- 
tibility for L = 24 and L = 32, respectively. The left inset 
in Figure |4] illustrates the scaling of the sample-to-sample 
variance of the sample average 



PSU = (i/Q)XX*- 



(20) 



Assuming that the width of these sample-to-sample fluc- 
tuations scales with the linear size L according to 



5([T*] av ) ~ L-^\ 



(21) 



we obtain, from the very good fitting for the sizes L = 
14 — 32 shown in the inset of Figure 01 the value v = 
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Fig. 5. FSS behavior of the ratio Ri w *] av (L) defined in equa- 
tion (|22[) for (a) the specific heat w — C and magnetic suscep- 
tibility w — x f° r h — 2 and (b) the magnetization w = M and 
the mean energy per spin w = (e) for h — 2.25. In panel (a) 
the solid lines show linear fitting extrapolations to the limit 
L — > oo. 



1.30(2). This estimate is in excellent agreement with the 
value v = 1.31(18) used in equation (fT5|) and estimated 
in our previous paper [36] from the scaling of the specific 
heat's pseudocritical temperature and consists further ev- 
idence for the breakdown of self-averaging in the RFIM. 

We turn now to study the behavior of the ratio R w 
defined above. In our formalism, we may define the vari- 
ance ratio R w in two distinct forms. One L-dependent 
ratio for the sample average of the individual maxima of 



w \w 



II',, 



„(£) 



V u , 



and one T-dependent for the average curve [w] 



R [w] av ( T ) - 



(22) 



(23) 



In Figure[5la) the behavior of the ratio R\ w *] av (£)> defined 
in equation (12"2"1) . where w* = \* — x{T — ^l,[x]* ) an d 
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Fig. 6. Low-temperature behavior of the ratio R[m] av (T) de- 
fined in equation (|23|l of the specific heat (w = C) and 
h = 1.9,2, and 2.1 in the main panel and the magnetic sus- 
ceptibility (w — x) an d h = 2 in the corresponding inset. The 
linear size dimension is L — 10. 



w * = C* = C{T = T LAC] ,J is illustrated for L = 16 - 32, 
as a function of the inverse linear size for disorder strength 
h = 2. The straight lines are extrapolations to L — > oc. 
In both cases a clear saturation to a non-zero limiting 
value is observed, i.e. Iimx-^oo R[x*] av (L) = 0.64(9) and 
liniL^oo R ^ c ,^ av (L) = .33(5), indicating violation of self- 
averaging |36I27I30|3"9"] . For the case of the susceptibility 
this violation is much stronger. Figure [5{b) presents the 
behavior of the ratio R\ w »] (L), where now w* = M* = 
M(T = T LM ,J and w* = (e>* = (e)(T = T i>[c] ,J as 
a function of the inverse linear size for sizes L = 4 — 32 
and disorder strength h = 2.25. It is clear that, while for 
the magnetization the ratio R\M*] av (L) approaches a non 
zero value in the limit L — > oo and therefore violation of 
self-averaging is expected, the same ratio for the mean en- 
ergy R[t e )*] av (L) starts with some very small values but 
finally approaches zero with increasing lattice size, indi- 
cating that in the infinite size limit the mean energy of 
the system is self-averaging. 

Finally, in Figure[5]we present the low-temperature be- 
havior of the ratio R\ w ] av (T), defined in equation (j23|) . for 
the specific heat w — C and the susceptibility w = \ (in- 
set) at h = 2 and for a lattice size L = 10. The ratio R[c] av 
for h = 1.9 and h — 2.1 is also shown in the main panel 
together with that of h — 2, calculated via the extrap- 
olation scheme described in reference |40j . Our intention 
is to identify the temperature variation of the non self- 
averaging property of the average specific heat and suscep- 
tibility. We observe that as we move from the disordered 
phase to the ordered phase this ratio increases, and gets its 
maximum value for the lowest temperatures shown. This 
indicates strong sample-to-sample fluctuations and non 
self-averaging behavior (which is more pronounced with 
increasing disorder strength), it becomes more evident as 
we enter the ordered phase, and is in agreement with pre- 
vious studies [25123118) . Since the sample-averaged pseud- 
ocritical temperature for L = 10 is of the order of [T*] oll ~ 



2.7, the violation of self-averaging goes well into the or- 
dered phase, and this should be compared with Figure 3 
of reference [39] • 

Let us comment here that, as stated above, possible 
breakdown of self-averaging may be traced back to situa- 
tions where the correlation length £ becomes of the order 
of the linear size of the system. The usual state of af- 
fairs is that the correlation length diverges just on the 
boundary between the ordered and disordered phase and 
consequently any break down of self-averaging may be ob- 
served only in the vicinity of the boundary [36 27 30 39 . 
However, as Dayan et al. [25] showed, the situation in the 
RFIM is quite different: the correlation length is of order 
of the linear size of the system everywhere in the ordered 
phase and not just at the boundary of the phase, indi- 
cating a violation of self-averaging in the ordered phase, 
as clearly manifested in our Figure |U Of course, self- 
averaging is expected to be restored at high-temperatures 
(above T c ) that £ -C L, because in this case, the system 
can be divided into independent regions of size £ and any 
measurement on the whole sample will be an average on 
these regions. This latter behavior is also depicted in our 
Figure^ where the ratio decreases rapidly and approaches 
zero for high temperatures for both the specific heat and 
the susceptibility. 

5 Conclusions 

In the present paper we investigated numerically some im- 
portant aspects of the critical properties of the d = 3 
RFIM with a bimodal random-field distribution in the 
strong disorder regime. In the first part of our study we 
focused on the scaling aspects of the model and presented 
evidence in favor of the so-called two-exponent scaling sce- 
nario. We also derived accurate values for the magnetic 
critical exponents j3, 7, and 7 of the model in good agree- 
ment with the existing estimates in the literature. 

In the second part, we investigated the self- averaging 
properties of the model in terms of probability distribu- 
tions of certain thermodynamic quantities and we illus- 
trated, using various finite-size measures and properly de- 
fined noise to signal ratios, the strong violation of self- 
averaging of the model in the ordered phase. Addition- 
ally, a FSS analysis of the width of the distribution of 
the sample-dependent pseudocritical temperatures of the 
magnetic susceptibility verified the theoretical expecta- 
tions of Aharony and Harris and provided an alterna- 
tive approach of extracting the value of the correlation 
length's exponent v. Overall, the data and analysis pre- 
sented in this paper indicate that the most complicated 
issue of breakdown of self-averaging in disordered systems 
deserves special attention and can be easily transformed 
into a useful tool that constitutes an alternative approach 
to criticality. 
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